% Generate Simulated stock returns from MJD model.
    % Output: logR(i,t) = simulated log return of market index under path i in year t. 
    % --> stored in file 'LogReturns.mat'

max_num_cont = 50;  % maximum number of future QRP contributions

%% Financial market framework

regime = "MJD";

r = 0.03;           % constant annual interest rate, compounded continuously.
infl = 0.02;        % annual rate of inflation.
MRP = 0.06;         % market risk premium (for investment in S&P 500, = 100% volatility).
delta = 0;          % dividend yield (assume =0 to make RILAs better comparable to TDFs).

if regime == "MJD"
    lambda = 0.20;      % likelihood of jump.
    mu_d_RN = 0.0735;   % mean of diffusion process (under risk-neutral measure).
    sigma_d = 0.14;     % volatility of diffusion process.
    mu_j = -0.25;       % mean of each jump.
    sigma_j = 0.10;     % volatility of magnitude of each jump.
elseif regime == "BS"
    lambda = 0;         % likelihood of jump.
    mu_d_RN = 0.03;     % mean of diffusion process (under risk-neutral measure).
    sigma_d = 0.18;     % volatility of diffusion process.
    mu_j = 0;           % mean of each jump.
    sigma_j = 0.01;     % volatility of magnitude of each jump.
else
    disp('Regime not specified correctly.')
    asdadsf
end

mu_d = mu_d_RN + MRP;   % mean of diffusion process (under real-world measure).

N_sim = 1e+6;   % number of simulated paths.
num_jumps = poissrnd(lambda,N_sim,max_num_cont);
logR = (mu_d - sigma_d^2/2)*1 + num_jumps*mu_j + sqrt( sigma_d^2 * 1 + num_jumps * sigma_j^2 ) .* randn(N_sim,max_num_cont);   % matrix of annual log returns, following MJD model.

save('LogReturns.mat','logR','r','infl')